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TRANSONIC AIRFOIL DESIGN 
USING CARTESIAN COORDINATES 

By Lei and A. Carlson 
Texas A&M University 

SUMMARY 

A numerical technique for designing transonic airfoils having a prescribed 
pressure distribution (the inverse problen) is presented. The method employs 
the basic features of Jameson's iterative solution for the potential equation, 
except that inverse boundary conditions and Cartesian coordinates are used. 

The method utilizes a direct-inverse approach that leads to a logical method 
for controlling trailing edge closure. Examples show the application of the 
method to design aft-cambered and other airfoils specifically for transonic 
flight. 

INTRODUCTION 

A numerical method for the design of transonic airfoils should not only 
be accurate but also be as simple as possible in concept and approach. It 
should use coordinate systems, input variables, and boundary condition treat- 
ments that can be easily understood by the user. In addition, it would be 
desirable if the method yielded the airfoil design shape for a given set of 
conditions with little or no iteration and used or computed nose and tail 
shapes that are aerodynamical ly and structurally reasonable. Finally, it 
should be able to handle both shocked and shockless flows and be suitable 
for not only complete design but also airfoil modification. 

Previous design methods and programs have either worked with the hodograph 
(1 21 (3l 

equations, ’ ' used direct optimization techniques, or tried the inverse 

( 4 - 7 ) 

approach. ' Of these, probably the best known are the hodograph methods of 

Nieuwland^^^ and of Bauer, Garabedian, and Korn.^^^ While these methods yield 

excellent accurate results, they are restricted to shock free solutions which 

may have adverse drag, pitching moment, and boundary layer separation charac- 

(8l 

ter i sties at off -design conditions.' Further, they are difficult from the 
user standpoint in that they involve complex mappings and transformations and 
require complicated initial input conditions. In Nieuwland's method, the user 



must a priori specify the trailing edge Mach number and flow angle, while in 
the Korn theory a set of logarithmic terms must be given involving singulari- 
ties inside the airfoil. As a consequence the hodograph methods may require 
a large number of trial solutions before an acceptable airfoil design is 
achieved, and they are not suitable for simple airfoil modification since 
there seems to be no straight forward way of changing only part of the solu- 
tion. 

Another design technique is to use a direct method (airfoil prescribed) 

to analyze the flow about a given airfoil and then, based upon this result, to 

modify the airfoil in an attempt to satisfy the design conditions. Usually 

this approach requires extensive experience on the part of the user and a large 

number of iterations. It can, however, yield good results. Recently, Hicks, 

( 3 ) 

Murman, and Vanderplaats have automated this procedure for non-lifting air- 
foils by combining a transonic, small perturbation analysis program with an 
optimization program based on the method of feasible directions. While very 
promising, the method tends to yield airfoils having multiple upper surface 
supersonic zones, which at off -design conditions may lead to multiple shocks 
and other adverse effects. In addition, it needs to be extended to lifting 
airfoils before its full usefulness can be ascertained. 

The last design formulation is the inverse method in which the airfoil 
surface pressures or velocities are specified and the airfoil shape subse- 
quently determined. Admittedly this approach requires knowledge of what would 
be a desirable pressure distribution, but this characteristic is probably un- 
derstood by the designer as well as, if not better than, any other. Further- 
more, the designer can a priori select a pressure distribution that will yield 
a given lift and moment j have a reasonable supersonic zone, behave desirably 
with, respect to boundary layer separation, and satisfy other transonic flow 
criteria. Of course, the resultant airfoil design shape may or may not be 
physically and structurally reasonable. 

The inverse approach was initially pursued by Erdos, Baronti, and ETzr 
weig^^^ and by Steger and Klineberg. Both groups used the transonic small 
perturbation equations and formulated the problem in a direct-inverse manner. 
That is, the leading edge shape was specified and that region solved directly, 
while on the remainder of the airfoil the pressure was prescribed and the shape 
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computed. Erdos et al used the perturbation potential as the basic unknown 
everywhere in the flowfield, but Steger and Klineberg used the perturbation po- 
tential only in the direct region. In the inverse zone, they posed the prob- 
lem in terms of the perturbation velocities. Since each approach solved the 
flow directly up to an arbitrary point, each could be used, if desired, to ana- 
lyze a problem in a completely direct mode. Now if the pressure distribution 
from such a direct computation were used as input for an inverse calculation, 
the original airfoil shape should be recovered. Yet, Erdos et al obtained in 
that case a discontinuity in body slope at the location of the shock wave, 
which they attributed to the singular nature of the point of intersection be- 
tween the shock and the surface boundary. Steger and Klineberg experienced 
similar phenomena but determined that they were numerical in origin and that 
they could be eliminated with careful formulation. Unfortunately, since the 
two methods used different dependent variables in the inverse zone, the self- 
consistency of an approach using only the perturbation potential was not estab- 
lished. 

This problem was subsequently investigated in detail^^^ and it was estab- 
lished that the discontinuity observed by Erdos was strictly numerical in ori- 
gin and that inverse results consistent with direct calculations could be ob- 
tained using the perturbation potential. In addition, it was demonstrated by 
using experimentally measured pressures as input that inverse transonic calcu- 
lations are capable of reflecting the consequences of viscous-inviscid inter- 
action on the airfoil displacement surface shape. Finally, it was determined 
that the inverse approach using the small perturbation equations tended to be- 
come inaccurate for thick and blunt-nosed airfoils. 

Quite obviously the latter problem can be eliminated by formulating the 
inverse problem with the complete potential flow equation. Unfortunately, this 
makes the problem more difficult in that unlike the small disturbance case the 
location of the airfoil surface boundary condition in the computational plane 
is unknown. Tranen^^^ applied the complete equations to the inverse problem 
by modifying the conformal mapping relaxation solution technique of Garabedian 
and Korn to use the pressure distribution boundary condition. Since this tech- 
nique maps the computational grid to match the airfoil surface, Tranen was 
forced to apply the boundary conditions at the surface of some assumed airfoil. 
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Only after the relaxation process had converged could a new shape be computed; 
and, usually, because of the difference between where the boundary conditions 
were and should have been applied, the resultant shape was not completely com- 
patible with the input pressure distribution. Tranen solved this problem by 
resorting to, for each case, a series of inverse-direct-inverse-direct, etc. 
calculations in which the pressure distribution was modified prior to each in- 
verse calculation. This was done in an attempt to achieve convergence, which 
occurred when the Cp distribution from a direct computation agreed with that 
used as input in the previous inverse run. 

Tranen 's work is significant because it successfully applies the complete 
equations to the inverse design problem and because it designs the entire air- 
foil. Since the method is entirely inverse and since in the inverse region 
the surface boundary condition specifies the derivative of (ji in the streamline 
direction, an initial value of <j) at the nose must be assumed. Unfortunately, 
the solution is sensitive to this value. Furthermore, to achieve convergence, 
Tranen' s approach requires large scale iteration involving an unspecified modi- 
fication of the input pressure distribution. 

The purpose of this report is to present and discuss a new numerical 
method suitable for the design and/or modification of subsonic and transonic 
airfoils. In order to achieve accuracy, the method utilizes the full inviscid 
potential flow equations; and in order to remain simple, it solves the problem 
in a stretched Cartesian grid system. The resultant working computer program 
has several unique features. First, it can be used in either the direct (anal- 
ysis) mode in which the airfoil shape is prescribed and the flowfield and sur- 
face pressures are determined, or in the direct-inverse (design) mode in which 
an initial nose shape is given along with the pressure distribution on the re- 
mainder of the airfoil and the flowfield and actual airfoil shape are computed. 
Second, it uses for the first time in a design program the rotated difference 
scheme introduced by South and Jameson, ^ which always has the correct 
zone of dependence in supersonic zones but does not require the coordinate sys- 
tan to be closely aligned to the flow direction. Third, unlike previous meth- 
ods, the present program determines the airfoil shape simultaneously with the 
flowfield relaxation solution. Thus, when the converged solution is achieved, 
an airfoil design compatible with the input pressure distribution is known. 
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and iteration is not required. 

This report primarily discusses the design mode of the program, presents 
the numerical features peculiar to the design problan, and shows how the 
method can be applied. A discussion of the analysis aspects of the method 
along with the numerical methods utilized, finite difference formulation, and 
numerical stability are presented in Reference 12. 

SYMBOLS 

a isentropic speed of sound 

a,b coordinate stretching constants 

c chord length 

Cp pressure coefficient, (p-P<„)/(^^U„^) 

C|_ lift coefficient 

Cf/|LE coefficient of moment about the leading edge 

f,g coordinate stretching functions 

M Mach number 

p pressure 

q velocity 

U,V velocity component in the x-, y- direction respectively 

x,y Cartesian coordinates 

a angle of attack 

Y ratio of specific heats 

r circulation 

e polar coordinate 

?,n computational coordinates 

p density 

$ potential function, Eq. (1) 

<j) perturbation potential, Eq. (2) 

Subscripts: 

«> freestream condition 

b body 

TE trailing edge 

i,j grid location 

3f 

C»n.x,y differentiation, i.e., fy = — 

^ oX 
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PROBLEM FORMULATION 


The exact equation for the potential function for two dimensional com- 
pressible flow can be written in Cartesian coordinates as 

(a^ - 4>x^)^xx “ 2$x^y^xy “ ®y^)'^yy ~ ^ 

where the subscripts denote partial differentiation. By defining a pertur- 
bation potential, ()>, such that 

$ = xqj:osa + yq„sinct + qj. 

where the velocity components are given by, 

U = $x = qjcosa + 

V = = qjsina + (})y) 

the governing equation for the perturbation potential becomes 

(a2 - u 2)4,^^ -2UV<j>xy + (a2 - v2)^yy = 0 (4) 

with 

[u2 + v2 - q.2] (5) 

The appropriate boundary condition at infinity is(13) 

<{) = ^Tan-'' (X^M^Tan(0 - a)) (6) 

where 9 is the polar angle, and r is the circulation, which is determined by 
the change in potential across the Kutta-Joukowski cut at the trailing edge 
of the airfoil , i .e. 

^ “ (‘*’y=0'*' " ‘**y=0~Hrail ing Edge 

As mentioned in the Introduction, in the design mode the shape of the 
nose region (typically 6-10% of the chord) is specified and the pressure is 
prescribed over the remainder of the airfoil. This procedure is used for sev- 
eral reasons. First, the nose region must be very accurately known in order 
to correctly fabricate an airfoil. Thus, by prescribing the nose shape, a 
possible major source of error is eliminated from the design procedure. Sec- 
ond, the boundary condition in the inverse region specifies the derivative of 
the perturbation potential in the tangential direction, and a starting value 


( 2 ) 

(3a) 

(3b) 
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must be known. With the present scheme, this value is determined by the di- 
rect solution in the nose region and need not be estimated or iterated for. 
Third, in some cases the designer may wish only to modify the aft portion of 
the airfoil. This can be done with the present scheme since the switch point 
from direct to inverse can be set by an input variable anywhere from about 6% 
chord to the trailing edge. Finally, and perhaps most importantly, specifi- 
cation of the nose shape gives the designer a physical entity whereby he can 
control the degree of closure of the trailing edge. This feature will be dis- 
cussed later. 

With this philosophy, the appropriate airfoil boundary condition in the 
direct region near the leading edge is 

b q^(sina+<|)y^) 

In the inverse region where the pressure is specified, the pressure coef- 
ficient expression can be solved for either <}>x or <|>y, since the tangential de- 
rivative involves both. If the inverse region included the leading edge re- 
gion, the logical choice would be to use (py near the nose, since the flow there 
is primarily in the y direction, and downstream. However, since in the 
present problem the nose is not normally part of the inverse region, the ap- 
propriate inverse boundary condition is 



Equation (9) is nonlinear in that (l)^, Ub^, and all involve values of (|) at 
various grid points. The treatment of this nonlinearity will be discussed in 
the next section. 

Equation (8) is applied on the assumed nose shape, and in the inverse re- 
gion Equation (9) is applied at the location of the current airfoil surfaces. 
The latter are determined by integration of the Equation (8) in the inverse 
region. 

To facilitate the application of the infinity boundary condition. Equation 
(6), the coordinate stretching represented on Figure 1 was selected. Here 
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the x-y plane is the physical plane and ?-n represents the computational plane, 
and each is subdivided into three regions. The stretching is symmetrical 
about the origin and is given by 

X = X4 + A2Tan[|-(c - C4)] + A3Tan[|-(c - (10) 

in region III and by 

X = ?(a + bg^) (11) 

in region II. The constants a and b are determined by the requirements 


X = X4 at C = C4 


(12a) 


and 


dx 

dC 



? = C4 


(12b) 


The constant A2 controls the grid spacing in the vicinity of X4, which is us- 
ually near the leading and trailing edges; while A3 determines the physical 
location of the grid adjacent to the grid edge. 

In the y-direction the stretching relationship is given by 


y = AiTan(|n) (13) 

where A-j controls the grid size near the airfoil. Here, only a one parameter 
stretching is used in order to facilitate the transformation from computational 
to the physical ordinates, which .is required each time the unknown airfoil 
surface is computed from Equation (8). 

Notice that these stretchings map the infinite x,y plane 

— 00 < < CO 

-00 <y<oo (14) 

into the finite computational plane 

- (1 + C4) 1 5 1 (1 + C4) 

- 1 < n 1 1 (15) 

where determines the amount of the computational plane confined to the vi- 
cinity of the airfoil. 

The governing equations formulated in terms of the ?-n variables are then 
written in finite difference form and solved iteratively using numerical 
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relaxation. The finite differences are formulated using a rotated scheme, 
first proposed by South and Jameson, ^ which always has the correct zone 
of dependence in supersonic regions but does not require the coordinate sys- 
tem to be closely aligned to the flow direction. Essentially, the finite 
difference expressions are arranged to exhibit the basic features of a local 
rotation to the streamline direction. 

With respect to the rotated difference scheme, it should be noted that 
the present approach is somewhat different from that used in Ref. (10) and 
(11). While still using rotation and viewing the relaxation process as a 
time-like procedure, time terms in the streamwise direction are not introduced 
implicitly as a consequence of the manner in which the difference expressions 
are formulated. Instead they are added explicitly and, as in Ref. (10) and 
(11), used as damping terms to control the stability and convergence of the 
relaxation process. By explicitly adding these time-like terms, no additional 
damping is required. Further, the amount of damping required can be easily 
determined by the user from the maximum local Mach number, which in the design 
problem would be known from the assumed surface pressure distribution. Also 
notice that the time-like terms correspond to the change in <t> between relaxa- 
tion cycles; and, thus, they approach zero as the solution converges. A de- 
tailed discussion of the finite difference formulation, the numerical scheme, 
and its stability is presented in Reference (12). 

SURFACE BOUNDARY CONDITIONS 
AND COMPUTATION OF SHAPE 

There are many ways to approximate the flow tangency. Equation (8), or 
surface pressure, Eq. (9), condition at the airfoil boundary. One approach, 
which is used here, is to generate dummy values of (f> at mesh points inside the 
boundary such that the usual difference equations can be solved at all points 
outside the boundary. The problem is to develop a scheme for providing and 
updating these dummy values by using the surface tangency or pressure condi- 
tion and neighboring values of <t> in the mesh, with adequate accuracy and with- 
out creating instability. 

To accomplish this for the direct region, first note that in the compu- 
tational coordinates Equation (8) becomes 
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( 16 ) 


^ _ Vb _ stng t 

Mx^b Ub cosa + fb<f> 5 ^ 


Here, Taylor series about the dummy point (see Figure 2) can be used 

to express and (|)^^. Thus, 

<l>n. = + (^K - ni.i 




n^i,j-l 

(17) 

>• • 

(18) 


When these are written in finite difference form using second order expres- 
sions for and and at least first order ones for and they become 
in the direct region (for the upper surface case) 




2a 


+ (n, - „j.,) ^ 


(19a) 


i - Vl.,1-1 ' *1-1.J-1 

?b 2 a^ 


(% 


.t) ( '*'i+1J ' ^’+l,j 


•l.,~ ^1-l.j 

2A5An 


+ <f> 


i~l >4-1 ) (19b) 


These expressions can then be substituted into Equation (16) and the result 
solved for a sufficiently accurate that is in terms of the neighboring 

potentials, body slope, and body position. 

In the inverse region the surface pressure boundary condition 


<(>xb = ■*'b'f'?b 


- -cos 







must be applied at the location given for the surface by the previous relaxa- 
tion sweep or last surface update. Thus, to start the inverse calculation an 
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initial profile is needed. While this profile should be reasonable and compat- 
ible with the shape selected for the nose region, it is not critical and the 
final airfoil may be quite different from the initial profile. 

As in the direct region, the Taylor series expression for Equation 
(18), is used in Equation (20); but in the inverse case considerable care must 
be used in selecting the finite difference forms for j_-j and <l>^Tii,j-l’ 

For example, since the resultant expression is to be solved for j_-[, one 
possibility is to use ^ and central differences about (i-l,j-l). While 
numerically stable, this approach has two disadvantages. First, since the 
dominant term in the expansion is (p^, central differencing causes to 

depend mostly upon Second, the pressure determining <(>i,j-l is 

Cp|^^ 1 "^ot Consequently, adjacent boundary points are uncoupled, 

and the boundary values usually have an undesirable oscillatory behavior. 

This can be mitigated slightly by using Cp^^. ^ and centering about (i-Js,j-l), 
but the results are still frequently oscillatory. 

In order to eliminate this oscillatory behavior, several other possible 
approaches were investigated. Among those found to be inaccurate and/or un- 


stable were first order forward schemes using Cp first order backward 
formulas using Cp|^^. , and schemes centered about (i,j-l) and using Cp,. The 
latter failure was expected since its algebraic form emphasizes (j)^ instead 
of (p^. Finally, it was determined that a scheme based on (i,j-l) and that 
uses Cp^., second order backward differences in c on both (p^ and (|>^^ and first 
order on n in <|)^^ frequently worked very well. 


However, the best form to use for in Equation (20) is (for the upper 


surface) 



+ (nb - rij_i) ( 



(21) 


where the pressure coefficient is Cp^.. This form does not introduce any nu- 
merical instabilities, builds into <l>i,j-i the upstream history of the airfoil. 


and uses C 


It leads to smooth boundary values and airfoil shapes with no 


apparent oscillations. 
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When Equation (21) is substituted into Equation (20) and the result solved 
for <()^- j_i» the final expression involves the square of the airfoil slope, 
Vb^/Ub^- Technically, this makes the equation nonlinear. However, since the 
surface slope is usually small in the inverse region, it has been found to be 
sufficiently accurate to hold (V/U)b constant at the value obtained for the 
last airfoil surface update. 

For those situations where the flow at point (i,j) is supersonic, the 
rotated finite difference scheme may require a value as well as 

Thus, in all cases a value of <t>i j _2 is determined by extrapolation as 

‘f’i j-2 ~ 

The above procedures for determining the values of the dummy mesh points 
inside the boundary is performed twice for the relaxation procedure of column 
i. First, either Eq. (16) or (20) are used with (ji values obtained from the 
previous relaxation sweep in order to obtain old values for the dummy points 
'^ijj-1 ‘f'i,j-2* Then, after the column i has been relaxed, they are used 

again with as many current values of as possible to obtain new values 
and <f|jj-2* manner, the dummy mesh points will have both old and new 

values just like regular mesh points, and they can be used directly in the 
finite difference formulas without special treatment. 

A similar procedure is used to satisfy the boundary conditions on the 
lower surface of the airfoil. 

Now in the inverse case the airfoil shape must also be computed. In 
other words, for a given (f> solution the differential equation given by Equa- 
tion (16) must be solved for the surface ordinates, y, as a function of x, 
with the initial conditions given by the surface slope and ordinate at the 
interface between the direct and inverse regions. In the present problem the 
location of the interface and the values of the initial conditions are known 
because the nose region is solved directly. If the problem were treated com- 
pletely inversely, the shape would depend upon the unknown value for (p at the 
leading edge. 

An easy method for solving Equation (16) is to use the Euler-Cauchy (or 
Runge-Kutta of order two) method, either directly or with iteration, since it 
only requires information at i and i+1 . This approach was tried, but it was 
found tc be insufficiently accurate. Consequently, the Runge-Kutta method of 
12 
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order four has been used; and while it requires for each step information at i 

i-Bs, and i+1, the improvement in accuracy is worth the added complexity. 

In the process of solving Equation (16) and (j)^^ must be evaluated by 

finite differences. While Equation (19a) is the obvious choice for (()„., 

'b 

either Equation (19b) or (21) could be used for (j> . Both have been used and 

both work well, but experience indicates that more accurate airfoil shape re- 
sults are obtained using the central difference scheme. Equation (19b). 

Finally, notice that in integrating Equation (16) the most difficult re- 
gion from an accuracy standpoint is near the leading edge where the surface 
slopes are very large. In the present approach this difficulty is essentially 
eliminated since the nose shape is prescribed and the nose region is solved 
directly. 


NUMERICAL STUDIES 

Self Consistency 

The present numerical method is capable of being used either directly or 
inversely. Thus, when pressure results from a direct calculation are used as 
an input for an inverse calculation, the original airfoil shape should be re- 
covered. However, as shown by References (4-6), this task can be difficult; 
the key to success is numerical consistency between the direct and inverse 
computations. Hence, this type of test is a check on not only the validity 
of the inverse approach, but also on the consistency of the difference schemes 
Several such tests have been conducted and some typical results are shown on 
Figures 3 and 4. Since the boundary condition in the inverse case is computed 
using Equation (20) combined with the backward difference formula, Eq. (21), 
the input Cp used in these tests (obtained from a direct calculation) was com- 
puted using Equation (21) for the term. In this manner, algebraic con- 
sistency between the two calculations is preserved. 

Figure 3 compares the actual surface slopes with those computed by the 
inverse technique for a nonlifting NACA 0012 at slightly supercritical con- 
ditions. In this case, the inverse design region extended from 6% chord to 
the trailing edge. As can be seen the slopes, and consequently the shape, are 
returned exactly. 
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A more severe test of self-consistency is a supercritical lifting case, 
and the results of such a test are presented on Figure 4. Here a strong shock, 
having a local upstream Mach number of about 1.30, exists on the upper surface 
of the NACA 0012 airfoil. Notice that the agreement between the direct and in- 
verse lift coefficients and between the actual and computed slopes is excel- 
lent, and that the differences between the upper and lower surface values are 
too small to be detected on the plot. For this case the computed surface or- 
dinates are everywhere within 0.33% (T/C) of the actual NACA 0012 ordinates, 

max 

and no problems exist near the shock or the trailing edge. 

Based upon these tests and results, it is believed that the present in- 
verse design scheme and program is self-consistent and reasonably accurate. 

Size of Inverse Region 

The present design approach is general in that the interface between the 
direct and inverse regions can be located anywhere between the leading and 
trailing edges. However, the nose region is very hard to compute accurately 
in the inverse mode if the surface slopes are steep; and the question arises— 
How close to the leading edge can the inverse region start? To answer this 
question a series of inverse calculations were performed for a six percent 
biconvex airfoil at of 0.9 and one degree angle of attack. The input Cp 
distribution, which contained shocks on both the upper and lower surfaces, 
was obtained from a direct computation, and in each of the tests the inverse 
zone was started at a different point. The inverse results were then compared 
with the actual surface slopes, given for the upper surface by the straight 
line = -0.24(x/c) + .12. 

In making these comparisons it should be noted that the -inverse boundary 
condition at point i, see Equations (20-21), uses information at points i, 
i-1, and i-2. Thus, if the inverse region starts at the third column from 
the leading edge, the inverse boundary condition uses only information aft of 
the leading edge. However, if the inverse region starts nearer the leading 
edge, information from in front of the airfoil will be used in the inverse 
boundary conditions, which may not be bad since the flow in front is subsonic 
and fully aware that there is an airfoil disturbing the flow. 
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The results of these tests are shown on Figure 5. The slopes computed 
by starting the inverse region at the third column from the leading edge are 
accurate, but those starting at the second are slightly in error on the slope 
and consequently also the shape. Finally, the slopes based upon the inverse 
region starting at the first column aft of the leading edge are in serious 
error. This trend has been verified by other tests with other airfoils. Con- 
sequently, it has been concluded that the inverse region should start no sooner 
than the third column from the leading edge, which is usually at about 5-10% 
chord. Hence, all the design results presented in this report have been com- 
puted with the inverse region starting no earlier than the third column aft 
of the leading edge. 

Sensitivity to Errors in Cp 

Since the pressure distribution is the primary input for the inverse mode, 
an important question is— How sensitive is the shape of the designed airfoil to 
errors in the input Cp distribution? This problem has been numerically in- 
vestigated, and it has been found that a 0.5% change in all the Cp input val- 
ues (ACp~0.003) leads to a 0.5% (T/C) change in the final airfoil ordinates. 
This increment appears small (A(y/c) ^ 0.0006), but it can be quite signifi- 
cant in the vicinity of the trailing edge where ordinates tend to be small. 
While this result is not at all unexpected, it does mean that care should be 
exercised when using the inverse design program to insure the accuracy of 
the input Cp distribution. 

COMPUTATIONAL PROCEDURE 

In designing an airfoil shape using the inverse scheme, the same basic 
numerical relaxation technique as described in Reference (12) is used. To 
start the problem an initial airfoil shape must be assumed, but this choice 
is not critical and the final airfoil shape may be considerably different. 

Next, since experience has indicated that the inverse scheme works best if 
the perturbation potentials have reasonable starting values, fifty relaxation 
cycles are performed in the direct mode for the initial airfoil shape on a 
very coarse grid. Then, the grid spacing is halved and the inverse procedure 
is initiated using the input pressure distribution as the boundary condition 
in the inverse region. 
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An important consideration in the inverse calculations is the frequency 
of computing or updating the new airfoil shape. Since in the early relaxa- 
tion cycles the pressure boundary condition may cause very large perturbations 
from the initial shape, it has been found best not to start the update pro- 
cedure until after the first fifty iterative cycles on the second grid, which 
are all in the inverse mode. After that, updating can occur after each re- 
laxation sweep, if desired, or at some specified interval, or after complete 
convergence has been achieved. The last is essentially the procedure used by 
Tranen.^^^ Since each update computation requires a finite amount of computer 
time, it would be inefficient to calculate a new airfoil shape after each re- 
laxation cycle. On the other hand, tests with the present scheme have indi- 
cated that updating only after convergence has been achieved usually yields 
airfoil shapes not completely compatible with the desired pressure distribu- 
tion. It appears that the optimum procedure is to perform the airfoil shape 
update calculation after every ten relaxation cycles, More frequent updating 
does not improve the accuracy of the final shape and only wastes computer 
time, while less frequent use sometimes leads to small errors in the final 
shape. After convergence on the second grid, the spacing is halved at least 
once more and the entire procedure repeated. 

Another important consideration in the inverse approach is the determina- 
tion of convergence. In direct calculations, it is usually sufficient to moni- 
tor both the maximum perturbation potential change between relaxation cycles 
and the number of relaxation cycles. The calculation is stopped when either 
the maximum perturbation potential change is less than some specified value or 
.the number of cycles exceeds a given number. However, in the inverse problem, 
another important criteria is the maximum change that occurs between update 
calculations in the surface ordinates of the airfoil being designed. It has 
been found that this value should also be monitored, and the computations con- 
tinued until it is consistently less than some specified value: (say Ay/c < 
0.0001). Frequently, satisfaction of this criterion yields a maximum pertur- 
bation potential change smaller than would be required for a direct calculation. 

All of the inverse results presented in this report have been obtained with 
this procedure and have used at least two grid halvings. Usually, the final 
grid contained 49 x 49 points, which yields 66 pressure points on the airfoil 
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surfaces. Typical computation times were six minutes on a CDC 6600 (including 
Varian plots) or about twenty minutes on an IBM 360/65. 

APPLICATION AND TYPICAL RESULTS 

When designing an airfoil to satisfy a specified pressure distribution, 
the final shape must have certain reasonable characteristics. Obviously, the 
upper and lower surface must not "criss-cross" near the trailing edge; and 
since the inverse approach yields the location of the boundary layer displace- 
ment surfaces, the trailing edge should have a finite thickness. Since the 
displacement thickness near the trailing edge is frequently on the order of 
0.4% of the chord, the computed trailing edge thickness should be on the order 
of 1.0% or more. Unfortunately, it is not at all obvious how the input pres- 
sure distribution should be modified to achieve such characteristics if they 
do not result. Besides, such a modification would change one of the desired 
design parameters, namely, the pressure distribution. 

With the present scheme, however, the pressure distribution can be kept 
the same, and the nose shape can be used by the designer to control the de- 
gree of tail closure. The procedure is demonstrated on Figure 6, which shows 
three airfoils, each designed with the same Cp distribution from 7% chord to 
the trailing edge. The input pressure distribution, which is shown as the 
solid line on Figure 7, is supposed to be "shockless" and contains the lower 
surface bucket typical of aft-cambered airfoils. 

Airfoil No. 4 has an NACA 0012 nose shape; but, as can be seen on Figure 
6, the resultant design has too thick a trailing edge. Consequently, a Korn- 
like nonsymmetrical nose shape (see Table I) having a smaller leading-edge 
radius was used. The resultant design, termed Airfoil No. 5 on Figure 6, is 
about two percent thinner and has a very reasonable trailing edge shape. Fin- 
ally, for Airfoil No. 6, the lower surface nose region ordinates were raised 
by 0.1% of the chord; and this change led to a slightly thinner airfoil with 
an even thinner trailing edge. 

Figure 8 shows another case in which the nose shape was used to control 
the trailing edge characteristics. The pressure distribution for these air- 
foils, which is the solid curve on Figure 9, was selected to have the same 
basic lift coefficient and lower surface pressure distribution as an NACA 0012 
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at the same flight conditions, but without the strong upper surface shock 
wave. In each case, the nose shapes and initial starting profiles were those 
associated with NACA OOXX airfoils. For the NACA 0010 nose shape. Airfoil No. 
100, the upper and lower surfaces of the designed airfoil criss-crossed; but, 
as can be seen on Figure 8 for Airfoils 110 and 115, as the leading edge rad- 
ius was increased, the maximum and trailing edge thicknesses also increased 
proportionately. In each case, the moment coefficient was low compared to 
similar aft-cambered airfoils, and the maximum upper surface ordinate was near 
the mid-chord. 

Based upon the results shown in Figures 6 and 8, it is believed that the 
nose shape can be used to control the trailing edge shape and that any desired 
thickness of the trailing edge displacement surfaces can be achieved by the 
designer by adjusting the nose shape. As indicated previously, this adjust- 
ment does not require changing the desired inverse Cp distribution. 

Now a severe test for a design program is whether or not an analysis or 
direct solution of the designed airfoil returns the design or inverse Cp dis- 
tribution. Figure 7 compares the inverse Cp used to design Airfoil No. 5 with 
that obtained from a direct solution (airfoil given) using the ordinates for 
No. 5, and Figure 9 makes a similar comparison for Airfoil No. 115. In these 
comparisons the input Cp is used with computed via backward differences 
(Equation (21)), while the direct Cp results are based upon more usual central 
differences. As a result, perfect agreement should not be expected. Never- 
theless, the agreement in both cases is excellent. On Figure 7, the lower sur- 
face result agrees perfectly; and the direct solution for the upper surface 
only disagrees slightly near 70% chord, exhibiting a weak shock (local Mach 
number of 1.07) at 84% chord. This appearance of a weak shock and the conse- 
quent variation in lift coefficient of about two percent frequently appears 

when aft-cambered airfoils that are designed to be shockless are analyzed di- 
( 2 ) 

rectly.' ' The comparisons for Airfoil No. 115 on Figure 9 show even better 
agreement and tend to verify the validity and accuracy of the present inverse 
design program. 

A final case is shown on Figure 10. An arbitrary pressure distribution 
(dashed line) having an upper surface Mach number plateau around 1.2 followed 
by a large jump to subsonic flow at 76% chord was selected for the inverse 
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input. No attempt was made to match the shock jump conditions at the dis- 
continuity. On the lower surface, the Cp was chosen to maintain subsonic 

flow; and a pressure bucket selected according to the Stratford separation 
(141 

criteria' ' was included to enhance lift. As indicated, the design program 

uses for a backward difference scheme, Eq. (21) with the C_ boundary con- 
b ^ 

dition, Eq. (20); and, thus, in regions of large gradients the output pres- 
sure distribution, which is computed using a central scheme for <j)p and should 

b 

be more accurate, will be different. 

Now in the course of the inverse solution, the trailing edge displacement 
surfaces that satisfied the input Cp distribution were not parallel; and, thus, 
the inviscid solution required a rear stagnation point behavior. The actual 
inverse pressure distribution, which was computed by central differences and 
is indicated on Figure 10 by the solid line, shows this behavior. (Of course, 
if the converged solution was used to compute Cp using backward difference 
for (j)r as in Eq. (21), the prescribed "inverse input" distribution would be 
obtained.) In addition, notice that the upper surface discontinuity has been 
smoothed. Examination of the flowfield results shows a smooth supersonic bub- 
ble on the upper surface and indicates that the decell eration, while rapid, 
is not due to a shock wave of any significant strength. Also shown on Fig- 
ure 10 are the results of a direct solution for the designed airfoil, which 
exhibit excellent agreement with the actual inverse pressure coefficient dis- 
tribution. Finally, the airfoil shown is the shape obtained after the boundary 


,( 2 ) 


for 


layer displacement thickness, calculated by the method of Nash-McDonald 
a Reynolds number of twenty million, has been subtracted. 

It is believed that the results of Figure 10 demonstrate that the present 
inverse technique and program can handle "arbitrary" pressure inputs. In ad- 
dition, it will yield results consistent with analysis technique. 


CONCLUDING REMARKS 

An inverse technique suitable for designing two-dimensional transonic 
airfoils having a specified pressure distribution has been developed. This 
method utilizes the full inviscid potential flow equation and exact boundary 
conditions and solves them via numerical relaxation. A rotated finite dif- 
ference scheme is employed in order to obtain the correct domain of dependence 
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in supersonic regions, and Jameson's time-like damping is included to insure 
numerical stability. Further, grid halving is used to achieve computational 
efficiency. 

A logical method for controlling trailing edge closure has been devel- 
oped. In addition, it has been demonstrated that it is not necessary to 
match the computational grid to the airfoil surface and that very accurate, 
numerically consistent, and physically correct results can be obtained in a 
Cartesian grid system. The technique has been successfully used to design 
aft-cambered and other airfoils specifically for transonic flight. 


Texas A&M University 
Aerospace Engineering Department 
College Station, Texas 
14 December 1974 
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Table I--Nose Shape Ordinates 
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Figure l--Flowfield Subdivision for 0 
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Figure 2--Relationship Between Airfoil and Grid (Upper Surface) 
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NOSE SHAPE NO. I T/C = O.IOI 


AIRFOIL NO. 5 



AIRFOIL NO. 6 
M^ = 0.79, «=0® 


Figure 6--The Use of Nose-Shape to Control Trailing Edge Closure, 
Example 1 
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Figure 7--Comparison of Inverse Cp Distribution with that Obtained 
by Analysis of Designed Airfoil, Example 1 
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Figure 8— The Use of Nose Shape to Control Trailing Edge Closure, 
Example 2 


31 












NASA-Langley, 1976 CR-2578 


33 



I 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
WASHINGTON. D.C. 20546 


OFFICIAL BUSINESS 
PENALTY FOP PRIVATE USE S300 


SPECIAL FOURTH-CLASS RATE 
BOOK 


POSTAGE AND FEES PAID 
NATIONAL AERONAUTICS AND 
• PACE ADMINISTRATION 
481 



J 

1 


I 


i 


POSTMASTER : 


If Undelfverable (Section 158 
Postal Manual) Do Not Return 


"The aeronautical and space activities of the United States shall be 
conducted so as to contribute . . . to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof." 

— National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS; Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 

Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. Also includes conference 
proceedings with either limited or unlimited 
distribution. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include final reports of major 
projects, monographs, data compilations, 
handbooks, sourcebooks, and special 
bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs, 
Technology Utilization Reports and 
Technology Surveys. 



Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION OFFICE 
NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 





